Vorticity Transport Equation

$$\frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + R_\nu ^{ - 1}{\nabla ^2}\omega$$

In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
from IPython import display
import math as mt
import matplotlib.animation as animation
from tqdm import tqdm
from numba.decorators import jit
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
def VTE(nx,ny,nt,cfl,mu,w,isav):
    global dt,dx,KX,KY,KX2,KY2,KXD,KYD
    
    dx=2*np.pi/nx;dy=2*np.pi/ny
    dt=0.01 #temporal time step
    
    ### define grids ###
    x=np.arange(nx)*dx;y=np.arange(ny)*dy
    X,Y=np.meshgrid(x,y)
    kx =2*np.pi/(dx*nx)*np.r_[nx/2,np.arange(1,nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/(dy*ny)*np.r_[ny/2,np.arange(1,ny/2),np.arange(-ny/2,0)]
    kxd=np.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)]   #de-aliased
    kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)]   #de-aliased
    kx2=kx**2; ky2=ky**2
    KX,KY  =np.meshgrid(kx ,ky )
    KX2,KY2=np.meshgrid(kx2,ky2)
    KXD,KYD=np.meshgrid(kxd,kyd)

    #w =+np.exp(-((X-np.pi+np.pi/10)**2+(Y-np.pi)**2)/0.1)\
    #   +np.exp(-((X-np.pi-np.pi/10)**2+(Y-np.pi)**2)/0.1)
    #   +np.exp(-((X-2*np.pi-np.pi/10)**2+(Y-2*np.pi-np.pi/5)**2)/0.1)
    wf=sf.fft2(w)

    whst=np.zeros((nt//isav,nx,ny))
    for it in tqdm(range(nt)):

        ### Cranck-Nicholson (only for the dissipation term)###
        #wf=1/(1/dt+0.5*mu*(KX2+KY2))*((1/dt+f(wf)+0.5*mu*(KX2+KY2))*wf)
        
        ### 4th-order Runge-Kutta ###
        g1=f(wf          ,1)
        g2=f(wf+0.5*dt*g1,0)
        g3=f(wf+0.5*dt*g2,0)
        g4=f(wf+    dt*g3,0)
        wf=wf+dt*(g1+2*g2+2*g3+g4)/6
                
        if(it%isav==0):
            w   =np.real(sf.ifft2(wf))
            whst[it//isav,:,:]=w         

    return locals()


def f(wf,ityp):
    phif=wf/(KX2+KY2)
    vxf = 1j*KY*phif; vx=np.real(sf.ifft2(vxf*KXD*KYD))
    vyf =-1j*KX*phif; vy=np.real(sf.ifft2(vyf*KXD*KYD))
    wxf = 1j*KX*wf  ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
    wyf = 1j*KY*wf  ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
    advw =vx*wx+vy*wy
    advwf=sf.fft2(advw)

    if(ityp==1): dt=cfl*dx/(max(np.amax(vx),np.amax(vy)))
    
    return -mu*(KX2+KY2)*wf-advwf
In [3]:
nx=256; ny=256; nt=12000; isav=100
mu=0.0002; cfl=0.4
w=np.random.randn(nx,ny)*5

dat=VTE(nx,ny,nt,cfl,mu,w,isav)
locals().update(dat)
100%|██████████| 12000/12000 [21:18<00:00,  9.39it/s]
In [4]:
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

def update_anim(it):    
    ax.clear()

    ax.imshow(whst[it,:,:],origin='lower',aspect='auto',cmap='jet')
    ax.set_title(r"$\omega$")
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim
Out[4]:

Vorticity Transport Equation in MHD

\begin{gathered} \frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + \left( {{\mathbf{b}} \cdot \nabla } \right)j + R_\nu ^{ - 1}{\nabla ^2}\omega \hfill \\ \frac{{\partial a}}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)a + R_\mu ^{ - 1}{\nabla ^2}a \hfill \\ \end{gathered}

In [5]:
def VTE_MHD(nx,ny,nt,cfl,mu,nu,w,j,isav):
    global dt,dx,KX,KY,KX2,KY2,KXD,KYD
    
    dx=2*np.pi/nx; dy=2*np.pi/ny
    dt=0.01 #temporal time step

    ### define grids ###
    x  =np.arange(nx)*dx
    y  =np.arange(ny)*dy
    kx =2*np.pi/(nx*dx)*np.r_[nx/2,np.arange(1,nx/2),np.arange(-nx/2,0)]
    ky =2*np.pi/(ny*dy)*np.r_[ny/2,np.arange(1,ny/2),np.arange(-ny/2,0)]
    kxd=np.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)]   #de-aliased
    kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)]   #de-aliased
    kx2=kx**2; ky2=ky**2
    X  ,Y  =np.meshgrid(x,y)
    KX ,KY =np.meshgrid(kx ,ky )
    KX2,KY2=np.meshgrid(kx2,ky2)
    KXD,KYD=np.meshgrid(kxd,kyd)
    
    #w = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
    #   -np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)\
    #   +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi-np.pi/5)**2)/0.4)
    #j = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
    #   +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)
    wf=sf.fft2(w)
    jf=sf.fft2(j)
    af=jf/(KX2+KY2)

    whst=np.zeros((nt//isav,nx,ny))
    ahst=np.zeros((nt//isav,nx,ny))
    jhst=np.zeros((nt//isav,nx,ny))
    for it in tqdm(range(nt)):

        ### 4th-order Runge-Kutta ###
        gw1,ga1=f(wf           ,af           ,1)
        gw2,ga2=f(wf+0.5*dt*gw1,af+0.5*dt*ga1,0)
        gw3,ga3=f(wf+0.5*dt*gw2,af+0.5*dt*ga2,0)
        gw4,ga4=f(wf+    dt*gw3,af+    dt*ga3,0)
        wf=wf+dt*(gw1+2*gw2+2*gw3+gw4)/6
        af=af+dt*(ga1+2*ga2+2*ga3+ga4)/6
        
        if(it%isav==0):
            w=np.real(sf.ifft2(wf))
            a=np.real(sf.ifft2(af))
            j=np.real(sf.ifft2(af*(KX2+KY2)))
            whst[it//isav,:,:]=w  
            ahst[it//isav,:,:]=a
            jhst[it//isav,:,:]=j
    
    return locals()

def f(wf,af,ityp):
    jf  =af*(KX2+KY2)
    phif=wf/(KX2+KY2)
    vxf = 1j*KY*phif; vx=np.real(sf.ifft2(vxf*KXD*KYD))
    vyf =-1j*KX*phif; vy=np.real(sf.ifft2(vyf*KXD*KYD))
    bxf = 1j*KY*af  ; bx=np.real(sf.ifft2(bxf*KXD*KYD))
    byf =-1j*KX*af  ; by=np.real(sf.ifft2(byf*KXD*KYD))
    wxf = 1j*KX*wf  ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
    wyf = 1j*KY*wf  ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
    axf = 1j*KX*af  ; ax=np.real(sf.ifft2(axf*KXD*KYD))
    ayf = 1j*KY*af  ; ay=np.real(sf.ifft2(ayf*KXD*KYD))
    jxf = 1j*KX*jf  ; jx=np.real(sf.ifft2(jxf*KXD*KYD))
    jyf = 1j*KY*jf  ; jy=np.real(sf.ifft2(jyf*KXD*KYD))
    advw =vx*wx+vy*wy
    advj =bx*jx+by*jy
    adva =vx*ax+vy*ay
    advwf=sf.fft2(advw)
    advjf=sf.fft2(advj)
    advaf=sf.fft2(adva)

    if(ityp==1): dt=cfl*dx/(max(np.amax(vx),np.amax(vy)))

    return -advwf+advjf-nu*(KX2+KY2)*wf,-advaf-mu*(KX2+KY2)*af
In [8]:
nx=256; ny=256; nt=12000; isav=100
mu=0.0002; nu=0.0002
cfl=0.4
w=np.random.randn(nx,ny)*5
j=np.random.randn(nx,ny)*5

data=VTE_MHD(nx,ny,nt,cfl,mu,nu,w,j,isav)
locals().update(data)
100%|██████████| 12000/12000 [41:19<00:00,  4.84it/s]
In [9]:
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(111)

def update_anim(it):    
    ax1.clear()#; ax2.clear()

    ax1.imshow(jhst[it,:,:],origin='lower',aspect='auto',cmap='jet')
    ax1.contour(ahst[it,:,:],80,linewidths=0.5,colors='k')
    ax1.set_title(r'$j(color\ map)+a(contour\ plot)$')#; ax2.set_title(r"$\omega$")
    ax1.set_xlabel(r'$x/(c/\omega_{pe})$')#; ax2.set_xlabel(r'$x/(c/\omega_{pe})$')
    ax1.set_ylabel(r'$y/(c/\omega_{pe})$')

anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)    
plt.close()
anim
Out[9]: